// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file
\brief	'C' Implementation of the Fast Epistasis (BT) Algorithm
*/
#include <math.h>

#include "fast_pepi_bt_assoc.h"
#include "pmaths.h"

double fast_pepi_assoc(int e1, BitString bs1, int e2, BitString bs2, struct sample *samples, int nSamples, int *summary_sig, int *summary_good, int *best_partner, double *best_score, double epi_alpha2, BOOL *is_valid) {
	double z = 0, or_aff, v_aff, or_unf, v_unf;
	register int i;
	long bitpos_bs1 = 0, bitpos_bs2 = 0;

	// Odds ratio test
	// make two 2x2 tables
	int a11 = 0, a12 = 0, a21 = 0, a22 = 0;
	int u11 = 0, u12 = 0, u21 = 0, u22 = 0;

	// Loop through all people in the sample.
	for (i = 0; i < nSamples; i++) {
		struct sample *person = samples[i].pperson;
		BOOL a1, a2, b1, b2;

		// Skipping person if phenotype missing, must increment the
		// bitpos counters to the next 2 'sample' records in the SNP bitstrings.
		if (person->missing) {
			bitpos_bs1 += 2;
			bitpos_bs2 += 2;
			continue;
		}

		a1 = BitStringgetbit(bs1, bitpos_bs1++);
		a2 = BitStringgetbit(bs1, bitpos_bs1++);
		b1 = BitStringgetbit(bs2, bitpos_bs2++);
		b2 = BitStringgetbit(bs2, bitpos_bs2++);

		// if affected 
		if (person->aff) {	      
			if (!b1) {
				//   ??x00
				if (!b2) {
					if (!a1) {
						// 00 x 00 
						if (!a2) a11 += 4; 
						else {
							// 01 x 00
							a11 += 2; 
							a21 += 2; 
						} 
					}
					else if (a2) a21+=4;  // 11 x 00
				}
				else  
				{
					//   ??x01
					if (!a1) {
						if (!a2) { // 00 x 01 
							a11 += 2; 
							a12 += 2; 
						}  
						else {  // 01x01
							a11++; 
							a21++; 
							a12++; 
							a22++; 
						} 
					}
					else if (a2) { // 11 x 01
						a21 += 2; 
						a22 += 2; 
					}  
				}
			}
			// ?? x 11 
			else if (b2) {
				if (!a1) {
					if (!a2) a12 += 4;    // 00 x 01 
					else { 
						// 01 x 01
						a12 += 2; 
						a22 += 2; 
					} 
				}
				else if (a2) a22+=4;  // 11 x 01
			}
		}
		// unaffected?
		else  {
			if (!b1) {
				//   ??x00
				if (!b2)  {
					if (!a1) {
						if (!a2) u11 += 4;   // 00 x 00 
						else { // 01 x 00
							u11 += 2; 
							u21 += 2; 
						}  
					} else if (a2) u21 += 4;  // 11 x 00
				} 
				//   ??x01
				else  {
					if (!a1) {
						if (!a2){ // 00 x 01 
							u11 += 2; 
							u12 += 2; 
						}  else { // 01x01
							u11++; 
							u21++; 
							u12++; 
							u22++; 
						} 
					} else if (a2) { // 11 x 01
						u21 += 2; 
						u22 += 2; 
					} 
				}
			}
			//  ?? x 11 
			else if (b2)  {
				if (!a1) {
					if (!a2) u12 += 4;  // 00 x 01 
					else { // 01 x 01
						u12 += 2; 
						u22 += 2; 
					}  
				} else if (a2) u22 += 4;   // 11 x 01
			}
		} 
	}

	// Calculate log(OR) and SEs
	or_aff = log((double)(a11 * a22)/(double)(a12 * a21));
	v_aff = (1/(double) a11) + (1/(double)a12) + (1/(double)a21) + (1/(double)a22);
	or_unf = log( (double)(u11 * u22)/ (double)(u12 * u21) );
	v_unf = 1/(double)u11 + 1/(double)u12 + 1/(double)u21 + 1/(double)u22;
	z = fabs((or_aff-or_unf) / sqrt(v_aff + v_unf));

	// Check this is a proper result
	if (realnum(z)) {  
		if (summary_good != NULL && summary_sig != NULL && best_score != NULL && best_partner != NULL) {
			// Count as a good result
			summary_good[e1]++;
			summary_good[e2]++;

			// Do we want to record this as part of the summary for the first set? 	  
			if (z >= epi_alpha2) {
				summary_sig[e1]++;
				summary_sig[e2]++;
			}

			// Is this result the best scrore yet for marker in set A?
			if (z > best_score[e1]) {
				best_score[e1] = z;
				best_partner[e1] = e2;
			}

			// The second marker might also be in set A 
			if (z > best_score[e2]) {
				best_score[e2] = z;
				best_partner[e2] = e1;
			}
		}

		*is_valid = TRUE;
	} else  {
		z = JUNK_Z; // Rubbish to assign -ve value to flag a rubbish result. 
		*is_valid = FALSE;
	}

	return z;
}